STUDY AREA

The main part of the Usa Basin (93500 km², Fig. 1) is located in the Komi Republic, except for some of the northern sections that extend up into the Nenets Autonomous Region. The Usa River discharges into the Pechora River on the west side of the catchment. The Ural Mountains with elevations ranging from around 300 to 1900 m bound the area in the east (Fig. 2a). To simplify the analysis, we restricted our scope only to the lowlands, as the processes controlling vegetation distribution in the mountains are to some extent different from those in the lowlands. Therefore, we excluded all areas above 250 m, after which our study area covers 79,500 km². The altitude of this lowland part of the basin ranges from 40 m a.s.l to mainly below 200 m a.s.l. Topography is generally very smooth, except for some river valleys that are up to several tens of meters deep.

There is a southwest- to northeast-oriented declining gradient in the mean annual temperature in the area. Meteorological data for the period 1961–1990 show that the mean annual temperatures varied from −2.5°C in Pechora (located immediately to the southwest of the Usa Basin) to −6.1°C in Vorkuta (located in the northeastern part of the Usa Basin). The climate is relatively continental; the mean January and July temperatures were −20.3°C and +16.1°C in Pechora, and −21.2°C and +13.0°C in Vorkuta, respectively. Mean annual precipitation in the region is around 550 mm, except for areas near the Urals where orographic precipitation results in higher values. About 75% of the Usa Basin is dominated by permafrost with various degrees of continuity, ranging from isolated patches in the south to continuous permafrost in the north. In the lowlands, permafrost temperature ranges from −4.5°C in the north to just below freezing point in the south (Fig. 2c, Oberman and Mazhitova, 2001, 2003).

Vegetation zones in the Usa Basin lowlands range from taiga in the south to forest-tundra and tundra in the north. The northern part is covered by tundra vegetation, the upland areas being occupied by dwarf shrub tundra vegetation with a well-developed lichen and/or moss layer. Peat plateau mires are common. Willow (Salix spp.) dominated, often paludified, vegetation occurs in depressions and river valleys. The central part of the Usa Basin consists of a mosaic of tundra and northern coniferous taiga forests, and the southern part belongs to the northern taiga forest zone (Kozubov et al., 1999). Large open mires are common in the lowlands of the taiga zone. The forest stands in lowland areas mainly consist of mixed forests dominated by spruce (Picea obovata Ledeb.). Downy birch (Betula pubescens Ehrh.) is the most frequent broadleaved tree. Forest line forests are also mainly spruce dominated, but typically some birches are found as well. Scots pine (Pinus sylvestris L.) is rare and is only found tens of kilometers to the south of the spruce treeline. A few larches (Larix sibirica Ledeb.) are encountered, especially on better drained soils.

Apart from the surroundings of a few villages, towns, and industrial areas, there has been almost no forest cutting. Also other human impacts on the vegetation, like construction of roads, have been minor, with the exception of vegetation changes around the industrial complex of Vorkuta (Virtanen et al., 2002) and some oilfields located especially north of Usinsk (Fig. 2a).

CLIMATE DATA

A regional climate model was developed for our study region by the Danish Meteorological Institute (more details in Christensen and Kuhry, 2000). We further downscaled this 14-yr data (1980–1993) from the original 16 km grid to a 1 km grid of the main climatic variables: monthly mean temperatures, degree-day sums (with threshold values of 0, +5, or +10°C), and minimum and maximum temperatures (the mean over the years 1980–1993 of the lowest and highest extreme yearly temperatures, see Fig. 2b). The degree-day sums, i.e. the sum of the daily mean temperatures over the threshold temperature, have been used to explain forest line and biome locations (Tchebakova et al., 1994; Tuhkanen, 1999).

Downscaling was made by calculating multiple regression models for climatic variables explained by latitude, longitude, and altitude. For degree-day sums, maximum temperature, and the mean monthly temperatures from April to October, the r² of the regression models varied from 0.94 to 0.99, and for the rest of the months and minimum temperature from 0.67 to 0.91. Precipitation of the HIRHAM model was downscaled to a 1 km grid by using kriging interpolation and adjusted precipitation data (explained in Van der Linden and Christensen, 2003).

The monthly mean temperatures are realistically simulated by the HIRHAM model (Christensen and Kuhry, 2000). However, the values of the HIRHAM model extreme minimum temperatures are relatively high and their range narrow (−40.5°C to −36.6°C). For example, in the Khoseda-Khard meteorological station (Fig. 2a), the measured mean yearly extreme minimum temperature for years the 1980–1993 was −45.7°C, but the model predicts for the same years and location −38.7°C. It is likely that even at 16 km resolution important aspects of the extreme Arctic boundary layer is not captured well by the model. This includes the simulation of low level inversions and the related cooling over snow covered surfaces (but see also Rinke et al., 1999; Dethloff et al., 2002). However, it seems that the spatial pattern that is most important for our forest line location models is relatively realistically reproduced when compared to available meteorological data, i.e., the extreme minimum temperature values recorded for the lowland forest line area from 1961 to 1990.

SOIL AND PERMAFROST DATA

Digital soil and permafrost data layers were compiled based on 1:200,000 topographic maps (Kuhry et al., 2002; Mazhitova et al., 2003; Oberman and Mazhitova, 2003). We downscaled the mean ground temperature (at a depth of zero annual temperature amplitude) to 1 km grid data from the permafrost map (Oberman and Mazhitova, 2003; Fig. 2c). Unfortunately the soil map was not feasible for our analyses as map polygons, by definition, usually included more than one soil class.

DIGITAL ELEVATION MODEL AND TWI

We created a Digital Elevation Model (DEM) with 100-m pixel resolution for the whole Usa Basin with ARC/INFO's TOPOGRID tool using the contours and hydrological layers in the Russian 1:200,000 digital maps as a source data (Fig. 2a). Using the DEM we calculated the Topographic Wetness Index, TWI (Fig. 2d; Beven and Kirkby, 1979, also called the Compound Topographic Index: Gessler et al., 2000). It is a quantification of the position of a site in the local landscape, and is defined as:

VEGETATION DATA

Vegetation classification based on Landsat TM 5 satellite images was produced with a resolution of 30 m for the entire Usa Basin. The classification was done by a supervised method using a mosaic of images, which was constructed from eight different images from five different dates. We used ground truth data and oblique aerial photographs taken in the summers of 1998, 1999, and 2000 in the classification and its accuracy assessment. According to accuracy test based on test points interpreted from photographs, the main vegetation types (forests, willow dominated stands and meadows, tundra heaths, mainly unvegetated areas, and water bodies) were separated to a relatively high accuracy: 84% of the test points were classified correctly. Details of the classification procedure and results are presented in Virtanen et al. (2004).

Pixels corresponding to >20 % tree crown cover in ground truth data were classified as forests. The proportion of the forests in different parts of the Usa Basin based on this classification is shown in Figure 2e. The original classification data included 21 vegetation type/land cover classes. For the modeling effort in this study, we regrouped these classes into three main groups: Forest (including classes Spruce forest, Spruce-fir forest, Pine forest, Mixed forest, Birch dominated stands, Forest cuttings), Tundra (Dwarf shrub moss tundra heath, Dwarf birch heath, Dwarf shrub lichen tundra, Sparse alpine tundra, Human impacted tundra) and Other (Willow stands, Meadows, Bog partly with few trees, Open bog, Wetland, Tundra with bare peat, Mainly bare land, Human infrastructures, Water bodies); see also Figure 3.

There are some limitations in determining the exact location of the forest line using Landsat TM images (Rees at al., 2002). Forest patches larger than a few hectares can be distinguished quite reliably, but as the grain size of the forest mosaic becomes smaller, single pixels (30 m × 30 m) usually cover more than one vegetation type. For example, isolated small spruce stands, occurring in river valleys and on small sandy hills even tens of kilometers north of the larger forest stands, are often smaller than the size of a pixel. Similarly, narrow strips of forests, found especially along rivers, usually become mixed up with the surrounding vegetation. An additional problem is caused by the spectral confusion between some forest types and areas of small-grained (often only some tens of square meters) peatland areas with admixtures of small ponds, birch, willows, and/or meadow vegetation. The problem of erroneous classification due to mixed pixels (the reflectance values of which are caused by several types of vegetation/landscape types within one pixel) is well known in satellite image classification studies (Campbell, 1996). Therefore, it is probable that most of the pixels classified as forests and located more than 10 to 20 km north from the forest line are misclassifications in our data (see Fig. 3). This conclusion is supported by our observations made on helicopter flights over the region, and hundreds of oblique aerial photos taken during these flights. To omit the potential confusion caused by those few misclassified forest pixels in the far north, we decided to use the 10% forest cover criterion in our modeling.

MODELING

For the analyses, we first calculated the proportion and mean patch size for each of the three main vegetation type classes (forests, tundra, other) for each grid cell using six different sample sizes: 1 × 1 km (n = 77102), 3 × 3 km (n = 8583), 6 × 6 km (n = 2128), 9 × 9 km (n = 937), 12 × 12 km (n = 517) and 25 × 25 km (n = 120). These were calculated using the program FRAGSTATS (McGarigal and Marks, 1995). Using the same sampling procedure, we calculated also the zonal mean values for the explanatory variables (climatic, soil, and topography variables). These data were imported to the SAS statistical program (SAS Institute, 1999), where logistic regression models for every scale, i.e. different grid size data, were constructed to predict the location of forests. Logistic regression allows one to predict a discrete outcome from a set of variables that may be either continuous, discrete or a mix. It is generally used to, for example, construct models connecting ecological phenomena to landscape-related variables (Franklin, 1995). To construct models with binary response, we classified grid cells with >10 % cover of forest pixels as forested, and the rest as nonforested. The used explanatory variables are continuous. The structure of the models is as follows:

FUTURE SCENARIOS

Developed logistic regression models allow us, by changing climatic input variables according to a chosen scenario, to predict the equilibrium location of the forest line under different temperature conditions. To illustrate the significance of even slight increase in temperatures, we calculated scenarios adding 1°C to July and minimum temperatures. Furthermore, we calculated for another scenario by adding 2.8°C to July and minimum temperatures; this corresponds to a predicted increase for our study region in the 2080s, based on the HadCM2S750 run of the second generation coupled ocean-atmosphere GCM model of the Hadley Centre (Mitchell et al., 1995; Johns et al., 1997).

ANALYSES LIMITATIONS DUE TO DATA RESOLUTION

Spatial resolution of various source data sets the lowest limit of scale that is possible to use in modeling (Fortin et al., 2000). The vegetation data has originally the highest spatial resolution of 30 × 30 m, while the climatic variables were statistically downscaled from 16 km resolution data to 1 km grid size. In the original ground temperature data the size of the smallest polygons was about 10 km². So, the downscaling of this data to 1 km grid cells generated some pseudo-resolution. It is also known that polygons in the permafrost map include small landforms of which ground temperatures lie outside the range indicated for the polygon.

LANDSCAPE STRUCTURE

According to our vegetation classification data, 25.8 % of the Usa Basin lowlands is covered by forests, 28.2% by tundra, and 46.0% by other cover types such as wetlands/peatlands (31.3%), willow-dominated areas or meadows (9.4%), water bodies (3.2 %), stony or sandy areas (1.6%), and mainly nonvegetated areas due to human impacts (0.5 %). The proportion of the forest-covered area gradually decreases from south to north. The width of the forest-tundra transition zone in the Usa Basin lowlands varies from about 100 to 140 km (Figs. 2e, 3). In addition to the variation due to the forest-tundra transition, especially extensive peatland areas but also willow dominated vegetation along the rivers and streams fragment both the tundra and forest landscapes (Fig. 3).

The fragmented landscape pattern can be seen also as a large variation in the mean patch size south from the forest-tundra transition zone in Figures 4a and b. On the other hand, there are no larger than 0.3 ha sized forest patches in areas where mean July temperatures are colder than about +13.3°C, or where ground temperature is below −1°C, corresponding to the area where more contiguous permafrost begins to prevail also in upland areas (Fig. 4b). The same type of pattern exists also in the forest area proportions. However, we decided to present the mean patch sizes here, as it illustrates one type of the landscape structure variables widely used in quantitative structure analysis, and which are potential input parameters for different models, such as the ones describing patch dynamics.

MODELS FOR THE FOREST LINE LOCATION

In logistic regression models, the most significant explanatory variable was the mean July temperature at all scales (Table 1). Although the degree-day sums were almost as good predictors, we only included the mean July temperature to the models, as all the variables describing summer temperatures are highly correlated. According to our model, the present forest line in the area roughly follows the isoline +13.9°C of the mean July temperature (Fig. 5a). The second significant variable, entered to every scale model, was mean ground temperature (Table 1). According to the models, forests do not exist in areas with ground temperature below 0.1°C (Fig. 5c). The third variable explaining significantly the location of the forest line at all scales except at 25 km was extreme minimum temperature (Table 1, Figs. 5b, 6). The effect of TWI was significant only at the 1- and 3-km scales (Table 1, Figs. 5d, 6). Other additional explanatory variables mentioned in Material and Methods section were not significant. The location of forests in our vegetation classification and the location predicted by the models at 1 × 1 km, 6 × 6 km, and 25 × 25 km resolution are presented in Figures 7a–c.

The relative significance of the explanatory variables varies between the scales (Table 1, Fig. 6). The mean July temperature and ground temperature seem to change nonlinearly and in opposite directions through different scales. The significance of the minimum temperature is higher at the more coarse scales, but TWI is significant only at finer scales. The degree of variation explained by the models rose from 44.9% at the finest resolution to 74.6% at the coarsest resolution model (Table 1). However, much spatial details are lost in predictions when coarser resolution models are used (Figs. 7a–c).

Based on our 1-km-scale model, presently 48% of the Usa Basin lowlands could be forested, but due to paludification, water bodies, etc., forests actually cover only 26% of the area. When mean July and minimum temperatures were increased by 1°C, potentially 77% of the Usa Basin lowlands could be forested. If we assume that areas presently occupied by land cover types other than forest or tundra will remain nonforested also in the future, the forest coverage would in the +1°C scenario be 43%, and forests could spread approximately 70 km northwards from their present location. When the climate scenario for the 2080s was applied to developed models, 99.8% of the Usa Basin lowlands would be climatically suitable for forest growth.

TOPOGRAPHIC AND CLIMATIC GRADIENTS BEHIND THE FOREST LINE LOCATION

The four different factors that significantly explained the location of the forest line in our models can be assumed to have different mechanisms connecting them to the presence of forests along the tundra-taiga transition zone. As our analyses are based only on correlative patterns, the actual ecological and physiological mechanisms behind the observed patterns discussed below are not more than tentative.

The importance of July mean temperature (or, some other variables describing summer temperature conditions) as a factor determining the location of the forest line can be well understood, as it indicates the amount of radiation energy available for growth. The mean July temperature isoline +13.9°C found in our study area for the forest line location is clearly higher than values of +10 to 12°C presented for the Siberian or North American regions (Timoney et al., 1992; Tuhkanen, 1999; MacDonald et al., 2000). Timoney et al. (1992) even indicated that the border of the continuous forest zone is found at +13°C in the central Canadian Arctic. Our threshold of 10 % forest cover in landscape is a stricter criterion than what has been used in some of the respective studies. With the same criterion, the coniferous forest line formed by Scots pine in northern Finland corresponds to +12.6°C mean July temperature isoline for the period 1961–1990 (Mikkola and Virtanen, unpublished).

One explanation for the exceptionally high temperature found in our study might be due to slight discrepancies between modeled climate data and meteorological observations in the area. The mean July temperature in 1980–1993 was +13.4°C in the Khoseda-Khard meteorological station, which locates at the forest line region in the middle of the Usa Basin (Fig. 2a). The HIRHAM model data used in our modeling gave the value of +13.8°C for the same location and years. It is also possible that the period used in developing climate models (1980–1993) was exceptionally warm when compared to earlier years. The mean July temperature in the Khoseda-Khard measured for the period of 1931–1960 was +12.5°C (data from 29 yr) and +13.4°C in 1961–1990. Thus, the differences due to above mentioned factors are not large enough to be the main explanation for the magnitude of differences found between our and earlier studies.

The location of forests is also explained by ground temperature. According to our models, forests are not found in areas where the mean ground temperature is below 0°C, which indicates the presence of permafrost. Permafrost has drastic effects on soil moisture conditions and the evidently related functioning and horizontal growth of roots, which again have effects on the water balance and mechanical stability of the trees. During our fieldwork, we did not observe spruce to grow on permafrost terrain (except for one case which could have been late seasonal ice). Possible rare exceptions might be forest stands on well-drained permafrost terrain with deep active layers. In North America and Siberia the situation is different, as there especially Picea mariana and Larix dahurica and also some other tree species grow on permafrost, and actually the thawing of ice-rich permafrost has destroyed tree stands (Camill, 1999; Osterkamp et al., 2000; Skre et al., 2002). Permafrost was largely, if not entirely, absent during the warm period of the Holocene in the Pechora lowlands (Oksanen et al., 2001; Väliranta et al., 2003), and it can be hypothesized that this has resulted in spruce ecotypes that have not adapted to grow on permafrost. This view is also partially supported by Kullman and Engelmark (1997) who suggested that ground frost conditions (and related risk of winter desiccation) are the most critical factors for spruce establishment and persistent growth in the treeline region of northern Fennoscandia.

While July temperatures explained the main south-north and altitude gradients, the extreme minimum temperature explained the pattern of forests locating further north (and in cooler summer temperatures) in the western parts of the Usa Basin compared to the eastern parts. The pattern indicated by minimum temperatures follows a climatic continentality gradient, which again reflects the level of variance or amplitude in climatic variables. The idea that it describes some large scale climate pattern related to the occurrence of forests is supported by the finding that the significance of minimum temperatures in the models is greater in coarser-scale models (Table 1, Fig. 6).

Edaphic features, like excessive soil moisture and the presence of extensive peatlands, also limit the occurrence of forest stands. The TWI derived from local topography is connected to soil moisture conditions, as elevated and often sandy areas have low TWI-values and are typically forested, but low and flat, often paludified areas having high TWI-values are often treeless or only sparsely stocked. As TWI was a significant explanatory variable only in finer-scale models (Table 1, Fig. 6), it indicates that some amount of variation in the location of forests can be explained only using spatially more detailed data, while small-scale variation is lost in models working at coarser scales.

PRESENT AND FUTURE FOREST LINE

Our simple scenario illustrated that even a slight increase in temperatures would allow extensive forest expansion into the present tundra. We observed in the field at some locations relatively young spruce individuals growing well in the former shrub tundra, located at least a few hundreds of meters from older and relatively large spruce stands. On the other hand, it seems that isolated small spruce stands typically located on small sandy hills, suggested relicts from earlier warmer periods of the Holocene, occurring even tens of kilometers north of the larger forest stands, have not produced new seedlings during the last decades (Lavrinenko and Lavrinenko, 1999).

It is probable that time lags exist in forest expansion due to different factors. Many studies have shown that at least under slightly warmer climatic conditions the growth and density of trees in stands existing at the forest line increase, but that the forests have not clearly expanded into tundra (Lavoie and Payette, 1994; MacDonald et al., 1998; Juntunen et al., 2002). Some studies, on the other hand, have reported relatively quick tree range expansion as a response to warming (Kullman, 2001), and fast responses are found especially in some paleoecological studies (Clark, 1998).

Plenty of potential factors can limit or at least delay the expansion of forests to new areas. Such factors in our study region might be quality of soil (nutrients, structure) and paludification (Sveinbjörnsson, 2000; Skre et al., 2002), permafrost and its changes (Oberman and Mazhitova, 2001; Van der Linden et al., 2003), fires (Rupp et al., 2000; Payette et al., 2001; Gromtsev, 2002), wind conditions (Sveinbjörnsson, 2000), reindeer (Oksanen et al., 1995), use of firewood (Hofgaard, 1999; Skre at al., 2002), and land-use changes due to extensive oil exploitation.

STATISTICAL MODELING AS A TOOL FOR UNDERSTANDING THE NATURE OF THE FOREST LINE

Models using variables that can be easily mapped or modeled and spatially extrapolated at the landscape or regional level are usable tools to understand and predict large biogeographic patterns like the location forest line. A multiscale approach, like the one presented in this study, can give more information about variables and their inter-relationships at various spatial scales (Levin, 1992; Franklin, 1995; Wu and Qi, 2000). Especially, our models give a more detailed insight into patterns occurring at finer spatial scales than the coarser scale global and large regions equilibrium vegetation models commonly used in climate change studies (Tchebakova et al., 1994; Sykes and Prentice, 1996; Kittel et al., 2000). As indicated by TWI in our model, topography related moisture conditions are important determinants of vegetation patterns mainly at smaller spatial scales, i.e. downwards from a few square kilometers. Therefore, one implication of our models is that this kind of variables should be included in modeling efforts aimed to work at these scales.

Although not assessed in this study, we stress that the extent of study area used in modeling is another important aspect of the scale as discussed, for example, in Fagan et al. (2003). Constructing the models with varying extents of the study area may have profound effects on the performance of models and the importance of different variables. Therefore, in further studies, it might be useful to test in more detail the effect of a varying extent of study area around the forest-tundra transition zone on the importance of different parameters.

The complexity of patterns and processes related to boundaries like forest lines, has called for the development of methods capable of quantifying the nature and the effects of change in these areas (Fagan et al., 2003). Especially, the dynamic modeling linked to statistically sound analysis and modeling of boundaries is one of the most important challenges in future studies (Rupp et al., 2000; Fagan et al. 2003). Still, the outputs of rather established modeling approaches, like the one described in this study, combined with landscape metrics, can provide valuable information on several fields of studies of processes related to forest line studies. These include, at least, monitoring the change of forest lines (Rees et al., 2002), planning of spatially effective sampling schemes and upscaling of ecological process models to the regional level (Kittel et al., 2000). When this kind of statistical models are coupled with GIS, they also enable spatially detailed future scenario analysis according to climate change scenarios (Franklin, 1995; Sykes and Prentice, 1996; this study).